Hardware-accelerated volume rendering methods for rectilinear grids include ray casting, texture slicing, shear-warp and shear-image rendering, and splatting, see, Roettger et al. “Smart hardware-accelerated volume rendering,” Eurographics/IEEE TCVG Symposium on Visualization, 2003, texture slicing, see, Rezk-Salama, et al., “Interactive volume on standard pc graphics hardware using multi-textures and multi-stage rasterization,” Proceedings of the ACM SIGGRAPH/EUROGRAPHICS workshop on Graphics hardware, pages 109-118, ACM Press, 2000, Engel, et al., “High-quality pre-integrated volume rendering using hardware-accelerated pixel shading,” Proceedings of the 2001 ACM SIGGRAPH/Eurographics Workshop on Graphics hardware, pages 9-16. ACM Press, 2001, Pfister et al., “The VolumePro real-time ray-casting system,” Proceedings of ACM SIGGRAPH 1999, pages 251-260. ACM Press, 1999, Wu, et al., “Shear-image order ray casting volume rendering,” ACM Symposium on Interactive 3D Graphics, pages 152-162, June, 2003, and Pfister, “The Visualization Handbook,” chapter Hardware-Accelerated Volume Rendering, Chris Johnson and Chuck Hansen (Editors), Academic Press, 2004.
Splatting reconstructs a continuous function by sampling a scalar field of values using 3D reconstruction kernels associated with each scalar value. For volume rendering, the continuous function is mapped to screen space as a superposition of pre-integrated 3D kernels, which are called 2D footprints, see, Westover, “Footprint evaluation for volume rendering,” Proceedings of ACM SIGGRAPH 1990, pages 367-376, August 1990. Splatting offers flexibility in terms of volume grids, including non-rectilinear, and mixing with point-sampled geometry, see Mao, “Splatting of non rectilinear volumes through stochastic resampling,” IEEE Transactions on Visualization and Computer Graphics, 2(2):156-170, 1996, and Zwicker, et al., “Ewa volume splatting,” IEEE Visualization 2001, pages 29-36, 2001.
Zwicker et al. describe a method for aliasing-free splatting. However, achieving interactive high quality EWA splatting is still difficult due to the computational complexity of EWA filtering.
Splatting benefits from the use of pre-integrated reconstruction kernels. Since Westover's original work, see above, and Westover, “Interactive volume rendering,” C. Upson, editor, Proceedings of the Chapel Hill Workshop on Volume Visualization, pages 9-16, May 1989, most volume splatting methods focus on improving image quality, including ray-driven perspective splatting, edge preservation, eliminating popping and blur, and image-aligned splatting, see, Mueller, et al., “Fast perspective volume rendering with splatting by utilizing a ray driven approach,” Proceedings of the 1996 IEEE Visualization Conference, pages 65-72, October 1996, Huang et al, “Edge preservation in volume rendering using splatting,” Proceedings of the 1998 IEEE symposium on Volume visualization, pages 63-69, 1998, Mueller et al., “Eliminating popping artifacts in sheet buffer-based splatting,” Proceedings of the 1998 IEEE Visualization Conference, pages 239-246, October 1998, Mueller et al., “Splatting without the blur,” Proceedings of the 1999 IEEE Visualization Conference, pages 363-370, October 1999, and Mueller et al., “High-quality splatting on rectilinear grids with efficient culling of occluded voxels,” IEEE Transactions on Visualization and Computer Graphics, 5(2):116-134, 1999.
With splatting, the volume data are interpreted according to 3D reconstruction kernels, one kernel for each particle, discrete sample point, or voxel. Each 3D reconstruction kernel has the effect of absorbing and emitting light. Integrals are predetermined separately across each 3D kernel, resulting in ‘footprint’ functions. Each footprint function ‘spreads’ the contribution of each point over nearby pixels in an output image. Typically, the span of the footprint is in the order of two to five pixels. The ‘smeared’ contributions of each voxel in the 3D volume are composited, i.e., accumulated, in a front-to-back or back-to-front order to produce the pixels of the output image.
Aliasing problems in volume splatting is described by Swan et al., “An anti-aliasing technique for splatting,” Proceedings of the 1997 IEEE Visualization Conference, pages 197-204, October, 1997, and Mueller, et al., “Splatting errors and antialiasing,” IEEE Transactions on Visualization and Computer Graphics, 4(2):178-191, April-June 1998. There, a distance-dependent stretch of the footprints was used to make the footprints act as low-pass filters.
Heckbert, in “Fundamentals in Texture Mapping and Image Warping” Master's Thesis, University of California at Berkeley, Department of Electrical Engineering and Computer Science, 1989, describes a rendering method that uses elliptical weighted average (EWA) filters to avoid aliasing of surface textures. However, that method only operates on regularly sampled 2D texture. In other words, that method cannot be used directly with irregular point samples or volume data sets.
Zwicker et al., extended Heckbert's method to represent and render texture functions on irregularly point-sampled surfaces, and to volume splatting, see Zwicker et al., “Surface splatting,” Proceedings of A CM SIGGRAPH 2001, pages 371-378, July 2001, Zwicker et al., “Ewa splatting,” IEEE Transactions on Visualization and Computer Graphics, 8(3):223-238, 2002, and Zwicker et al., “Ewa volume splatting,” IEEE Visualization 2001, pages 29-36, 2001.
Point-based models have been successfully rendered on a graphics processing unit (GPU), see, Rusinkiewicz et al., “Qsplat: A multiresolution point rendering system for large meshes,” Proceedings of ACM SIGGRAPH 2000, pages 343-352, July 2000, Botsch et al., “High-quality point-based rendering on modern GPUs,” Proceedings of the 2003 Pacific Graphics Conference, pages 335-343, 2003, and Guennebaud et al., “Efficient screen space approach for hardware accelerated surfel rendering,” Proceedings of the 2003 Vision, Modeling and Visualization Conference, pages 1-10, November 2003.
Ren, et al., in “Object-space ewa surface splatting: A hardware accelerated approach to high quality point rendering,” Proceedings of Eurographics 2002, pages 461-470, September 2002, described an object space formulation for EWA surface splats and an efficient implementation of the formulation on graphics hardware. For each point in object space, quadrilaterals, which were texture-mapped with a Gaussian texture, were deformed to result in a correct screen-space EWA splat after projection.
Other prior art method includes opacity-based culling, fast splat rasterization, hierarchical splatting, object and image space coherence, shell splatting, 3D adjacency data structure, and post-convolved splatting, see, Mueller et al., “High-quality splatting on rectilinear grids with efficient culling of occluded voxels,” IEEE Transactions on Visualization and Computer Graphics, 5(2):116-134, 1999, Huang et al., “Fastsplats: Optimized splatting on rectilinear grids,” Proceedings of the 2000 IEEE Visualization Conference, pages 219-227, October 2000, Laur et al., “Hierarchical splatting: A progressive refinement algorithm for volume rendering,” Proceedings of ACM SIGGRAPH 1991, pages 285-288, August 1991, Ihm et al., “On enhancing the speed of splatting with indexing,” Proceedings of the 1995 IEEE Visualization Conference, pages 69-76, 1995, Botha et al., “Shellsplatting: Interactive rendering of anisotropic volumes,” Proceedings of 2003 Joint Eurographics—IEEE TCVG Symposium on Visualization, May 2003, Orchard et al., “Accelerated splatting using a 3d adjacency data structure,” Proceedings of Graphics Interface 2001, pages 191-200, September 2001, and Neophytou et al., “Post-convolved splatting,” Proceedings of the symposium on Data visualisation 2003, pages 223-230, Eurographics Association, 2003.
Lippert et al., in “Fast wavelet based volume rendering by accumulation of transparent texture maps,” Proceedings of Eurographics 1995, pages 431-443, September 1995, described a splatting method that directly uses a wavelet representation of the volume data. A hierarchical and frequency sensitive splatting method is based on wavelet transformations and pre-computed splat primitives, which accomplished view-dependent and transfer function-dependent splatting, Welsh, et al., in “A frequency-sensitive point hierarchy for images and volumes,” Proceedings of the 2003 IEEE Visualization Conference, Seattle, USA, October 2003. None of those methods have been implemented completely on a GPU.
Some GPU-accelerated splatting methods use texture mapping hardware for the projection and scan-conversion of footprints, see Crawfis et al., “Texture splats for 3d scalar and vector field visualization,” Proceedings of the 1993 IEEE Visualization Conference, pages 261-266, 1993, and Botha, et al., above.
EWA Volume Splatting
Volume splatting interprets volume data as a set of particles that are absorbing and emitting light. To render the data, line integrals are precomputed separately across each particle, resulting in 2D footprint functions or splats in the image plane. The splats are composited back-to-front to determine the output image. Particles are represented by 3D reconstruction kernels A common choice is 3D elliptical Gaussian kernels. The notation Gv(t−p) represents an elliptical Gaussian kernel centered at a 3D point p with a 3×3 variance matrix V:
                                          G            V                    ⁡                      (                          t              -              p                        )                          =                              1                                                            (                                      2                    ⁢                    π                                    )                                                  3                  /                  2                                            ⁢                                                                  V                                                                    1                  /                  2                                                              ⁢                                    ⅇ                                                -                                      1                    2                                                  ⁢                                                      (                                          t                      -                      p                                        )                                    T                                ⁢                                                      V                                          -                      1                                                        ⁡                                      (                                          t                      -                      p                                        )                                                                        .                                              (        1        )            
In theory, Gaussian kernels have infinite support. In practice, the kernels are truncated to a given cutoff radius r, i.e., the kernels are evaluated only in a range(t−p)TV−1(t−p)≦r2,  (2)where usually 1≦r≦3. Furthermore, the choice of Gaussians as 3D kernels guarantees a closed-form footprint function after integration along viewing rays.
However, during perspective rendering, the splatting process usually results in aliasing artifacts due to a changing of sampling rate along the viewing rays while rendering. EWA volume splatting minimizes the artifacts by convolving the footprint function with a 2D low-pass filter, which yields an aliasing-free footprint function called the EWA volume resampling filter.
Zwicker, et al., in “Ewa volume splatting,” IEEE Visualization 2001, pages 29-36, 2001, described a closed-form representation of the EWA volume resampling filter that is based on the following two assumptions. First, the low-pass filter takes the form of a 2D Gaussian reconstruction kernel function. Second, a non-linear perspective transformation, which maps reconstruction kernels to image space, is linearly approximated using a Jacobian matrix.
To summarize the derivation of the EWA volume resampling filter, a rotational part of the viewing transformation that maps object space to camera space coordinates is given by a 3×3 matrix W. Denote camera space coordinates by u=(u0, u1, u2). Normalized, the origin of the camera space u=0 is at a center of projection and an image plane is a plane at a distance u2=1. Camera space coordinates of a voxel k are given by uk. Image space coordinates are denoted by x, and the image space position of voxel k is xk. Further, the Jacobian matrix of the perspective projection at a point uk in camera space to image space is a 3×3 matrix Juk:
                              J                      u            k                          =                              (                                                                                1                                          u                                              k                        2                                                                                                              0                                                                      -                                                                  u                                                  k                          0                                                                                            u                                                  k                          2                                                2                                                                                                                                          0                                                                      1                                          u                                              k                        2                                                                                                                                  -                                                                  u                                                  k                          1                                                                                            u                                                  k                          2                                                2                                                                                                                                                                                    u                                              k                        0                                                                                                                                  (                                                                              u                                                          k                              0                                                                                ,                                                      u                                                          k                              1                                                                                ,                                                      u                                                          k                              2                                                                                                      )                                                                                                                                                                              u                                              k                        1                                                                                                                                  (                                                                              u                                                          k                              0                                                                                ,                                                      u                                                          k                              1                                                                                ,                                                      u                                                          k                              2                                                                                                      )                                                                                                                                                                              u                                              k                        2                                                                                                                                  (                                                                              u                                                          k                              0                                                                                ,                                                      u                                                          k                              1                                                                                ,                                                      u                                                          k                              2                                                                                                      )                                                                                                                                        )                    .                                    (        3        )            
Given the 3×3 variance matrix Vk″ of a reconstruction kernel k in object space, its transformation through camera space to image space is Vk=JukWVk″WTJukT.
The EWA volume resampling filter ρk(x) is obtained by integrating the reconstruction kernel in image space along viewing rays. The reconstruction kernel is convolved with the Gaussian low-pass filter, which, as shown by Zwicker et al, above, yields the 2D footprint function
                                                        ρ              k                        ⁡                          (              x              )                                =                                    1                              2                ⁢                π                ⁢                                                                        J                                          u                      k                                                              -                      1                                                                                        ⁢                                                                        W                                          -                      1                                                                                        ⁢                                                                                                                                                  V                          ^                                                k                                            +                                              V                        h                                                                                                                      1                    2                                                                        ⁢                          ⅇ                                                -                                      1                    2                                                  ⁢                                                      (                                          x                      -                                              x                        k                                                              )                                    T                                ⁢                                                      M                    k                                    ⁡                                      (                                          x                      -                                              x                        k                                                              )                                                                                      ,                                  ⁢        where                            (        4        )                                          M          k                =                                            (                                                                    V                    ^                                    k                                +                                  V                  h                                            )                                      -              1                                .                                    (        5        )            
There, Vh is a 2×2 variance matrix of the Gaussian low-pass filter, which is usually chosen to be the identity matrix. The 2×2 variance matrix {circumflex over (V)}k is obtained by skipping the third row and column of Vk1. A matrix with a hat symbol, ‘^’, denotes the result of skipping its third row and column.
The problem with the prior art filtering is that after a particular filter is selected, the same filter is used for the entire volume data set, without considering how the volume voxels were generated, and without considering how the data set is sampled along the viewing rays during the generation of the output image.
It is desired to solve this problem.